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A TOPOLOGY-INDEPENDENT  SHAPE  MODELING  SCHEM 
By 

KSSSSEb. 


The  resulting  equation  of  motion  is  solved  by  employing  entropy-satisfying  upwind 
finite  difference  schemes.  We  also  introduce  a new  algorithm  for  rapid  advancement 
of  the  front  using  what  we  call  a narrow-band  vpdalion  scheme.  This  leads  to  sig- 
nificant improvement  in  the  time  complexity  of  the  shape  recovery  procedure  in  2D. 
An  added  advantage  of  our  modeling  scheme  is  that  it  can  easily  be  extended  to 
any  number  of  space  dimensions.  The  efficacy  of  the  scheme  is  demonstrated  with 
numerical  experiments  on  low  contrast  medical  images.  We  also  demonstrate  the 
recovery  of  3D  shapes. 


CHAPTER  1 
INTRODUCTION 


An  important  goal  of  computational  vision  is  to  recover  the  shapes  of  objects  in 
2D  and  3D  from  various  types  of  visual  data.  One  way  to  achieve  this  goal  is  via 
modd-based  techniques.  Broadly  speaking,  these  techniques,  as  the  name  suggests, 
involve  the  use  of  a model  whose  boundary  representation  is  matched  to  the  image 
to  recover  the  object  of  interest.  These  models  can  be  either  rigid  or  nonrigid.  In 
the  former  case,  we  have  what  are  popularly  known  as  correlation-based  template 
matching  techniques,  while  the  later  involves  a dynamic  model  fitting  process  to  the 
data.  Shape  models  serve  the  purpose  of  an  intermediate  stage  in  object  recognition 
tasks,  since  they  provide  a more  stable  and  useful  description  than  the  original  raw 
data.  In  this  dissertation,  we  present  a new  dynamic  approach  to  shape  modeling 
which  retains  the  most  attractive  features  of  existing  methods,  and  overcomes  their 
prominent  limitations.  For  the  rest  of  this  dissertation  we  will  use  the  word  shape 


model  for  a he 


■ representation  of  surfaces. 


Shape  reconstruction  typically  precedes  the  symbol 
The  shape  models  must  aid  in  recovery  of  detailed  at 
ing  only  the  weakest  of  the  possible  assumptions  about  the  observed  shape.  Several 
variational  shape  reconstruction  methods  have  been  proposed  and  there  is  abundant 
literature  on  the  same  |5,  37,  42, 6,  50, 22].  Generalised  spline  models  with  continuity 
constraints  arc  well  suited  for  fulfilling  the  goals  of  shape  recovery  (see  [6,  7,  41)). 
Generalized  splines  are  the  key  ingredient  of  the  dynamic  shape  modeling  paradigm 
introducer]  by  Terzopoulos  el  u/.,  [45],  Incorporating  dynamics  into  shape  modeling 


enables  ihc  crcalion  of  realistic  animation  in  computer  graphics  applications  and  for 
tracking  moving  objects  in  computer  vision.  Following  the  advent  of  the  dynamic 
shape  modeling  paradigm,  there  was  a flurry  of  research  activity  in  the  area,  with 
numerous  application  specific  modifications  to  the  modeling  primitives,  and  exter- 
nal forces  derived  from  data  constraints  [17,  46,  47,  49,  51,  10,  12].  However,  the 
aforementioned  schemes  for  shape  modeling  have  two  serious  limitations  - the  de- 
pendence of  the  final  surface  shape  on  the  initial  guess  made  to  start  the  numerical 
reconstruction  procedure,  and  a strong  assumption  on  the  object’s  topology.  The 
first  of  these  deficiencies  stems  from  the  fact  that  the  nonconvex  energy  functionals 
used  in  the  variational  formulations  have  multiple  local  minima.  As  a consequence 
of  this  feature,  the  numerical  procedures,  lor  convergence  to  a satisfactory  solution, 
require  an  initial  guess  which  is  “reasonably"  close  to  the  desired  shape.  Existing 
shape  representation  schemes  have  an  additional  shortcoming  in  that  they  lack  the 
ability  to  dynamically  sense  the  topological  changes  that  may  occur  in  the  shape  of 
interest  during  the  shape  recovery  process.  In  this  dissertation,  we  present  a mod- 
eling scheme  that  makes  no  assumption  about  the  object's  topology,  and  leads  to  a 
numerical  algorithm  whose  convergence  to  the  desired  shape  is  completely  indepen- 
dent of  the  shape  initialisation.  Our  method  can  also  recover  shapes  whose  topology 
changes  over  time,  c.g.,  the  cell  boundary  in  cell  division  [35]. 

The  framework  of  energy  minimisation  has  also  been  used  successfully  in  the  past 
for  extracting  salient  image  contours  - edges  and  lines.  Kass  ef  of.  [17]  used  energy- 
miniillixing  “snakes"  that  exhibit  an  affinity  toward  image  features  such  as  edges 
points  and  edge  segments,  while  maintaining  piecewise  smoothness.  The  weights  of 
the  smoothness  and  image  terms  in  the  energy  functional  can  be  adjusted  for  different 
kinds  or  behavior.  Snakes,  also  referred  to  as  active  contour  models,  arc  restricted 
examples  of  the  more  general  techniques  of  matching  deformable  surface  models  to 
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Figure  LI.  Teat  bed  for  our  topology-independent  shape  modeling  scheme, 
image  data  by  means  of  energy  minimization  (45].  The  scheme  seeks  to  design  energy 


to  high-level  processes.  In  the  absence  of  a well-developed  high-level  mechanism  to 
make  a choice  among  these  solutions,  an  interactive  approach  is  used  to  explore  the 
alternatives.  By  adding  suitable  energy  terms  to  the  minimization,  the  user  pushes 
the  model  out  of  a local  minimum  toward  the  desired  solution.  In  the  problem  area 
of  automatic  segmentation  of  noisy  images,  snakes  perform  poorly  unless  they  are 
placed  close  to  the  preferred  shapes.  In  a move  to  make  the  final  result  relatively 
insensitive  to  the  initial  conditions,  Cohen  [9]  defines  an  inflation  force  on  the  active 
contour.  This  force  makes  the  model  behave  like  an  inflating  balloon.  The  contour 
model  with  the  above  change  will  be  stopped  by  a strong  edge  and  will  simply  pass 


Although  the  inflation  force  prevents  the  curve  from  getting  “trapped"  by  iso- 
lated spurious  edges,  the  active  contour  model  cannot  segment  complex  shapes  with 
significant  protrusions  like  the  one  shown  in  figure  1.1(b).  Moreover,  despite  a good 


initialization,  tin;  active  contour  model,  due  to  its  arc-longlll  and  curvature  mini- 
mization properties,  cannot  be  forced  to  extrude  through  any  significant  protrusions 
that  a shape  may  posses.  One  possible  solution  to  this  problem  is  to  embed  the 
snake  model,  which  is  an  instance  of  a ID  thin-plate-membrane-splinc,  in  an  adap- 
tive environment  wherein  the  material  parameters  controlling  the  relative  strengths 
of  elasticity  and  rigidity  are  adapted  (see  [36)).  The  merits  of  such  an  approach 
are  suspect  since  it  is  not  always  possible  to  derive  criteria  upon  which  to  base  the 

protrusions  in  complex  structures.  Terzopoulos  [45]  has  shown  that  multiple  instances 
of  deformable  models  are  required  to  handle  shapes  with  several  distinct  parts.  This 
can  be  very  cumbersome,  for  it  involves  excessive  user  interaction  and  presumes  that 
the  shape  has  already  been  segmented  into  its  constituent  parts.  Instead,  we  propose 
a dynamic  shape  modeling  method  that  will  start  with  a single  instance  of  the  model 
and  will  “sprout"  branches  during  the  evolutionary  process.  Once  the  shape  has  been 
segmented  from  the  image,  its  constituent  part  structure  can  be  inferred  using  the 
algorithm  described  in  [18]. 


be  known  before  the  shape  recovery  can  commence.  However,  it  is  not  always  possi- 
ble to  specify  the  topology  of  an  object  prior  to  its  recovery.  As  a result,  most  shape 
recovery  schemes  make  strong  assumptions  about  the  object  topology.  Unknown 
topology  is  an  important  concern  in  object  tracking  and  motion  detection  appli- 
cations where  the  positions  of  object  boundaries  are  tracked  in  an  image  sequence 
through  time.  During  their  evolution,  these  closed  contours  may  change  connectiv- 
ity and  split,  thereby  undergoing  a topological  transformation.  A heuristic  criterion 
for  splitting  and  merging  of  curves  in  2D  which  is  based  on  monitoring  deforma- 
tion energies  of  points  on  the  elastic  curve  has  been  discussed  in  [35].  Returning  to 


static  domain,  more  recently  molecular  dynamics  has  been  used  to  model  surfaces 
of  arbitrary  topology  [40]-  Smoothness  and  continuity  constraints  are  imposed  by 
subjecting  a particle  system  to  interaction  potentials  which  locally  prefer  planar  or 
spherical  arrangement.  Particles  can  be  added  and  deleted  dynamically  to  enlarge 
and  trim  the  surface,  respectively,  while  the  system  dynamics  strive  continually  to 
organise  the  particles  into  smooth  shapes.  The  result  is  a versatile  method  with 
applications  in  surface  fitting  to  sparse  data  and  3D  medical  image  segmentation. 

priori  assumption  about  the  object’s  topology  can  be  made.  A single  instance  of  our 
model,  when  presented  with  an  image  having  more  than  one  shape  of  interest  (see 
figure  1.1(c)),  has  the  abifity  to  split  freely  to  represent  each  shape  (25,  26].  We  show 
that  by  using  our  approach,  it  is  also  possible  to  extract  the  bounding  contours  of 
shapes  mlh  holes  in  a seamless  fashion  (see  figure  6.9).  Finally,  we  apply  our  method 
to  3D  images  and  recover  the  surface  shapes  of  3D  objects  (see  figure  (6.12)). 

1.1  Qveiview 

In  this  section  we  briefly  outline  our  scheme  to  model  complex  shapes.  This 
method  is  inspired  by  ideas  first  introduced  in  Osher  and  Sclhian  (30,  39)  to  model 
propagating  fronts  with  curvature-dependent  speeds.  Two  such  examples  are  flame 
propagation  and  crystal  growth,  in  which  the  speed  of  the  moving  interface  normal 
to  itself  depends  on  transport  terms  modified  by  the  local  curvature.  The  challenge 

front  which  will  accurately  approximate  these  highly  unstable  physical  phenomena. 
Sothian  [38]  has  shown  that  direct  parameterization  of  the  moving  front  may  be 
unstable  since  it  relies  on  local  properties  of  the  solution.  In  contrast,  a method 
which  preserves  the  global  properties  of  the  motion  is  sought.  Osher  and  Sethian 


[30,  39]  achieve  this  by  viewing  the  propagating  surface  as  a specific  level  set  of  a 
higher-dimensional  function.  The  equation  of  motion  for  this  function  is  reminiscent 
of  an  initial  valued  Hamilton-Jacobi  equation  with  a parabolic  right-hand  side  and  is 
closely  related  to  a viscous  hyperbolic  conservation  law. 

In  our  work,  we  adopt  these  level  set  techniques  to  the  problem  of  shape  recovery. 
To  isolate  a shape  from  its  background,  we  first  consider  a closed,  nonintersecting, 
initial  hypersurface  placed  inside  (or  outside)  it.  Following  the  above  level  set  ap- 
proach, this  hypersurface  is  then  made  to  flow  along  its  gradient  field  with  a speed 
F(K),  where  K is  the  curvature  of  the  hypersurface.  As  in  [30],  we  adopt  a global 
approach  and  view  the  (Af  — 1)  dimensional  moving  surface  as  a level  eel  of  a lime- 
dependent  funclion  0 of  Af  space  dimensions.  The  equations  of  motion  written  for  this 
higher  dimensional  function  are  then  amenable  to  stable  enlropp-salisfying  numerical 
schemes  designed  to  approximate  hyperbolic  conservation  laws.  Topological  changes 
can  be  handled  naturally  in  this  approach,  since  a particular  level  set  {0  = 0}  of 
the  function  (fr  need  not  be  simply  connected.  However,  there  are  two  problems  that 
need  to  be  surmounted  before  we  can  use  this  design  for  shape  recovery.  First,  it 
is  required  that  we  stop  the  hypersurface  in  the  neighborhood  of  the  desired  shape. 
We  do  this  by  synthesizing  an  negative  speed  function  from  the  image.  Secondly,  we 
have  to  construct  an  eilension  of  this  speed  function  to  other  level  sets  (0  = C } 
in  order  to  give  a consistent  meaning  to  the  image-based  speed  term  at  all  points  in 
the  image  (see  figure  3.1).  In  the  following  sections  we  outline  a possible  solution  to 
these  problems. 

We  note  that  this  work  on  interface  motion  and  hyperbolic  conservation  laws  as 
discussed  in  [30,  38,  39],  has  been  applied  in  the  area  of  computer  vision  for  shape 
characterization  by  Kimia  el  at.  [IS,  19],  who  unify  many  diverse  aspects  of  shape  by 
defining  a continuum  of  shapes  (roaction/dilfusion  space),  which  places  shapes  within 


a neighborhood  of  oilier  similar  shapes.  This  leads  to  a hierarchical  description  of 
a shape  which  is  suitable  for  its  recognition.  The  key  distinguishing  feature  of  our 
work  from  that  of  Kimia  et  of.,  is  that  they  assume  the  boundary  of  the  object  shape 
to  be  known,  while  we  reconstruct  it  from  noisy  data.  In  other  words,  they  show 
that  by  evolving  a known  shape  boundary , explicit  clues  can  be  derived  towards  the 
goal  of  developing  a hierarchical  shape  description.  In  contrast,  we  start  with  an 
arbitrary  function  if  and  recover  complex  shapes  by  propagating  it  along  its  gradient 
field.  Shape  characterisation  can  be  done  once  the  object  shape  is  extracted. 

An  important  contribution  of  this  dissertation  over  our  previous  work  in  [26] 
is  the  design  and  implementation  of  a faster  numerical  technique  for  solving  the 
governing  equation  in  2D  Our  new  algorithm  exploits  the  fact  that  the  front,  which 
is  a particular  level  set  [if  = 0}  of  a higher  dimensional  function,  can  be  advanced 
by  updating  the  function  if  at  a small  set  of  points.  This  scheme  is  an  alternative  to 
updating  the  if  function  at  all  the  points  in  the  computational  domain.  The  set  of 
points  at  which  the  update  procedure  is  applied  belong  to  a narrow  band  lying  on 
either  side  of  the  level  set  [0  = 0].  Since  the  narroio-band  update  strategy  involves 
only  a fraction  of  the  total  number  of  points,  a significant  saving  in  time  is  realised 
making  our  method  a very  attractive  alternative  to  other  shape  recovery  schemes.  A 
complete  discussion  of  the  narrow-band  techniques  for  interface  propagation  may  bo 
found  in  [1].  We  also  show  that  the  surface  shapes  of  objects  can  be  recovered  from 


In  summary,  we  present  a novel  scheme  for  shape  modeling  which  can  be  used  in 
both  computer  vision  and  computer  graphics  applications.  Given  the  reconstructed 
shape,  our  approach  can  also  be  extended  to  decipher  the  constituent  part  structure 


for  high-level  processing.  The  remainder  of  this  dissertation  is  organized  as  follows: 
chapter  2 introduces  the  curvature-dcpendcnt  front  propagation  problem  and  estab- 
lishes a link  between  Hamilton-Jacobi  equations  and  a hyperbolic  conservation  law. 
In  chapter  3 we  explain  our  level  set  algorithm  for  shape  recovery  and  in  chapter  4 
we  outline  the  details  of  our  narrow-band  approach.  Finally,  in  chapter  5 we  present 
some  experimental  results  of  applying  our  method  to  some  noisy  and  low  contrast 
medical  images  in  2D  as  well  as  3D.  We  close  with  a discussion  of  advantages  of  our 


ENERGY  METHODS  FOR  SHAPE  MODELING 


In  this  chapter,  wc  review  the  some  of  the  literature  on  existing  shape  modeling 
methods  in  computer  vision.  Specifically,  wc  provide  the  details  of  two  “energy- 
based"  methods  for  shape  modeling  and  reconstruction.  The  first  is  the  symmetry- 
seeking  deformable  model  [45)  and  the  second  is  a shape  recovery  scheme  in  2D  called 
active  contour  model  or  “snake"  [17], 

Existing  shape  modeling  methods  roughly  fall  under  two  categories:  passive  and 
active  models.  Models  of  shape  such  as  generalized  cylinders,  introduced  by  Binford 
[4),  and  lumped-parametor  family  of  shapes  such  as  superquadric  models  [3,  32]  are 
purely  geometric,  hence  passive.  Generalized  cylinders  arc  used  to  model  elongated 
shapes  with  axial  symmetry,  while  the  superquadric  shape  models  are  well  suited 
for  object  recognition  tasks  because  one  can  express  them  compactly  using  a small 

Lagraugian  dynamics  and  react  to  external  forces  in  a way  that  is  similar  to  real 
elastic  objects.  Generalized  spline  models  with  continuity  constraints  [7,  10,  12,  41, 
44,  49,  51]  are  prime  examples  of  the  active  shape  modeling  paradigm. 


The  deformable  model  is  a physically  based  modeling  framework  for  shape  and 
motion  reconstruction  of  free  form  flexible  objects.  In  this  framework  the  objects 
are  modeled  as  elastically  deformable  bodies  subject  to  continuum  mechanical  laws. 
Constraints  are  expressed  as  forces  which  deform  the  elastic  models  and  propel  them 


through  potentially  complicated  motions  such  that  they  satisfy  the  available  con- 
straints over  time.  There  are  two  types  of  forces,  intrinsic  and  extrinsic. 


while  the  extrinsic  constraints  reflect  observations  about  the  environment  that  can 
be  extracted  from  sensory  data.  The  model  is  embedded  in  a force  field  which  en- 
codes the  available  sensory  information.  The  ambient  force  field  will  then  mold  the 
deformable  model  to  make  their  3D  shape  consistent  with  the  observed  data.  A 
plethora  of  external  force  fields  can  be  synthesized  from  the  data  using  a wide  variety 
of  constraints.  Contrary  to  the  conventional  models  of  3D  shape,  the  deformable 
models  are  active  in  that  they  react  to  extrinsic  forces  as  one  would  expect  real 
elastic  objects  to  react  to  physical  forces.  This  is  because  the  deformable  models 
arc  governed  by  the  principles  of  elasticity  theory  as  expressed  through  Lagrangian 
dynamics.  Specifically,  the  continuum  mechanical  equation 
fl’v  9v  SEIv)  ,,  , 

"aF+'a+-ir  = f'v>  <21> 

governs  the  nonrigid  motion  of  the  deformable  model  in  response  to 
force  f(v),  where  Ji  is  the  mass  density  function  of  the  deformable  body  and  7 Is 


We  begin  this  subsection  by  providing  an  informal  description  of  the  deformable 
model.  Consider  a deformable  sheet  made  of  elastic  material  (a  mcinbrane-thin-platc 


material  is  passed  down  the  length  of  the  tube.  At  regularly  spaced  points  along  the 
spine,  it  is  coupled  to  the  tube  with  radially  projecting  forces  so  as  to  maintain  the 
spine  in  approximate  axial  position  within  the  tube.  Additional  forces  arc  included 
that  coerce  the  tube  into  a quasi-symmclric  shape  around  the  spine.  Finally,  extra 


control  is  provided  over  the  shape  by  introducing  expansion/conlraclion  forces  ra- 
diating from  the  spine.  The  rigidities  of  the  spine  and  the  tube  arc  independently 
controllable,  and  their  natural  rest  metrics  and  curvatures  can  either  be  prescribed 
or  modified  dynamically.  If  the  circumferential  metric  of  the  tube  is  set  to  zero,  for 
instance,  the  tube  will  tend  to  contract  around  the  spine,  unless  the  other  forces  pre- 
vail. The  model  will  shorten  or  lengthen  as  the  longitudinal  metrics  of  the  tube  and 
the  spine  are  modified.  In  short,  a variety  of  interesting  behaviors  can  be  obtained 
by  adjusting  the  control  variables  designed  into  the  model. 

The  deformable  tube  and  the  spine  can  be  described  by  geometric  mappings 
from  material  coordinate  domains  into  Euclidean  three-space  JP,  The  mapping  is 
expressed  by  position  vectors  in  space  whose  component  functions  are  time  varying. 
The  spine  is  a deformable  space  curve  defined  by  mapping  a univariate  material 
coordinate  domain  s e [0, 1]  into  JP 

v : [0,1]  x (0,  oo]  JP  such  that  v(s,  t)  = (X(s,  1),  /(s,  (),  Z(s,  I))  (2.2) 

The  tube  is  constructed  from  a deformable  sheet  that  is  defined  by  mapping  a bivari- 
ate material  coordinate  domain  (u,  o)  € [0,  l]2  into  Up 

v : [0,1]*  x [0,co]  -r  JP  nuch  that  v(u,u,t)  = (X(u,u,t),  K(u,u,t),  Z(u,u,t))  (2.3) 

The  tubular  shape  is  created  by  specifying  boundary  conditions  on  two  opposite 
edges  of  the  sheet  which  effectively  glue  these  edges  together.  The  edges  u = 0 and 
u = I are  “seamed”  together,  letting  a span  the  length  of  the  tube.  The  required 


v(l,o,t), 


(2.4) 

(2.5) 


(2.6) 


The  strain  energy  £s  associated  with  the  spine  function  v(s,t)  is  given  by 

£S(v)=/o 

where  the  vertical  bars  denote  the  Euclidean  norms  and  the  quantities  id,, is 
both  functions  of  (s,t).  <",(a,  i)  determines  the  local  tension  along  the  spine  and 
“if3.1)  determines  Us  local  rigidity.  The  weighting  functions  in,, to,  control  the 
elastic  properties  of  the  spine.  We  can  regulate  the  equilibrium  shape  of  the  spine  by 
suitably  defining  the  weighting  functions. 

The  strain  energy  £T  associated  with  the  tube  function  v(u,»,()  is  given  by 

£T(v) = j:  i ■» 1 ^ i’ + ^ i % i*  ■ +■» 1 0 11 + 

2,011  m 

The  weighting  functions  ui,o  and  too,  define  the  natural  metric  of  the  tube  along 
each  parameter  curve.  These  functions  locally  control  the  tension  of  the  deformable 
sheet  (that  constitutes  the  tube)  along  each  material  coordinate  curve.  Analogous 
expressions  for  uiu(u,o),ui„(tl,v)  and  too,(u,o)  determining  the  natural  curvature 
of  the  sheet  can  be  written  down.  These  functions  locally  control  the  rigidity  of  the 
sheet.  The  continuity  of  the  deformable  spine  can  be  controlled  via  the  functions 
in,  and  in,.  Position  discontinuity  is  introduced  by  setting  hi, (s,t)  = to,(s,i)  = 
0 and  tangent  discontinuity  by  setting  in,(s,l)  = 0.  Similarly,  we  can  introduce 
discontinuities  in  position  as  well  as  tangents  for  the  deformable  tube. 

The  spine  and  the  tube  functions  are  coupled  by  first  identifying  that  n = s, 
bringing  into  correspondence  the  spine  coordinate  with  the  coordinate  along  the 
length  of  the  tube.  The  mapping  function  of  the  spine  vs(e,  ()  is  then  distinguished 
from  that  of  the  lube  vT(n,u,t)  with  superscripts  S and  T. 


Hill 


(2-1)  wsociitc-d  with  tl 


The  force . 


The  snake  is  a model  of  a deformable  contour  composed  of  abstract  clastic  ma- 
terials. Two  types  of  materials  determine  the  snake  configuration:  string  and  rod. 
The  former  makes  the  snake  resistant  to  stretching  and  the  latter  makes  it  resis- 
tant to  bending.  This  deformable  curve  is  then  embedded  in  a potential  functional 
synthesised  from  a variety  of  image  constraints.  Depending  on  the  nature  of  the 
external  potential  field,  the  snake  will  extract  different  image  feature.  Unlike  most 
other  techniques  for  finding  salient  image  features,  the  snake  model  is  active.  It  is 
always  minimising  its  energy  functional  and  therefore  exhibits  dynamic  behavior. 

The  basic  snake  model  is  a amtroltal-contmuihj  spline  [41)  under  the  influence  of 
image  forces  and  external  constraint  forces.  The  internal  spline  forces  serve  to  impose 
a piecewise  smoothness  constraint  and  the  external  constraint  forces  are  responsible 
for  pulling  the  snake  near  the  desired  local  minimum. 

Consider  a deformable  curve  v(s,t)  with  parameters  s and  i (time),  defined  on 
given  intervals  11  — [0..1],  and  T,  respectively.  Let  this  deformable  curve  be  a function 
of  two  coordinate  variables  x and  y with  the  same  parameterization: 

v(s,t)  = (x(s,t), »(»,!)) : s e fl,(  6 T.  (2.14) 

The  potential  energy  functional  of  the  snake  £„„it(v)  is  defined  as 

ftnah(v)  = jf  (ft,.  + E„,)ds,  (2.15) 

where  Era  is  the  potential  associated  with  external  constraint  forces.  denotes 
the  internal  potential  functional  of  the  snake  and  is  given  by: 

Eia(v)  = jf  u>,(s)  | v,  I1  +tcj(s)  | v„  I’  ds.  (2.J6) 

The  first-order  term  to,(s)  | v.  |2  is  responsible  for  C“  continuity  in  the  snake  and  the 
second-order  term  tic,(.s)  | v„  |2  enforces  Cl  continuity.  The  weight  u>i(s)  regulates 


the  tension  of  the  snake,  whereas  10,(5)  regulates  its  rigidity.  Position  and  tangent 
discontinuities  may  be  introduced  along  a snake  by  setting  these  weights  to  aero. 

Shape  recovery  in  2D  can  be  accomplished  by  defining  a simple  energy  functional 

£=<(v)  = -Ja  | 70,  */(*,»)  p,  (2.17) 
then  the  snake  is  attracted  by  the  local  minima  of  the  potential,  which  corresponds 
to  the  local  maxima  of  the  image  gradient.  Here,  G„  * /(x,y)  denotes  the  image 
function  convolved  with  a Gaussian  smoothing  filter  whose  characteristic  width  is  a. 
Potential  functionals  modeling  other  external  constraints  defined  by  the  user  may  be 

Given  the  total  potential  energy  of  the  snake  (see  equation  (2.  IS))  for  a spe- 
cific initial  position,  the  equilibrium  configuration  can  be  readied  by  minimizing  it. 
Specifically,  if  the  function  v-  represents  a local  minimum  for  £,„t„  it  satisfies  the 
the  following  Euler-Lagrange  equation: 

r>;  - £(u>.(s)v.-)  + ^(Uri(s)v;.)  = -VE«,(v‘).  (2.18) 
In  this  formulation,  each  term  appears  as  a force  applied  on  the  curve.  Given  the 
value  of  v(s,l  = 0),  the  initial  snake  configuration,  the  above  equation  can  be  solved 
numerically  by  adopting  direct  or  iterative  finite  difference  schemes  [17,  20). 


CHAPTER  3 

FRONT  PROPAGATION  PROBLEM 


One  numerical  approach  lo  this  problem  is  to  take  the  above  Lagrangian  descrip- 
tion of  the  problem,  produce  equations  of  motion  for  the  position  vector  x(»,l),  and 
then  discretize  the  parameterization  with  a set  0r  discrete  marker  particles  laying 
on  the  moving  front.  These  discrete  markers  are  updated  in  time  by  approximat- 
ing the  spatial  derivatives  in  the  equations  of  motion,  and  advancing  their  positions. 
However,  there  arc  several  problems  with  this  approach,  as  discussed  in  Sethian  [38], 
First,  small  errors  in  the  computed  particle  positions  are  tremendously  amplified  by 
the  curvature  term,  and  calculations  arc  prone  to  instability  unless  an  extremely  small 
time  step  is  employed.  Second,  in  the  absence  of  a smoothing  curvature  (viscous) 
term,  singularities  develop  in  the  propagating  front,  and  an  entropy  condition  must 
be  observed  to  extract  the  correct  weak  solution.  Third,  topological  changes  are  dif- 
ficult to  manage  as  the  evolving  interface  breaks  and  merges.  And  fourth,  significant 
bookkeeping  problems  occur  in  the  extension  of  this  technique  to  three  dimensions. 

•3,2  Level  Set  Formulation 

As  an  alternative,  the  central  idea  in  the  level  set  approach  of  Osher  and  Sethian 
[30]  is  to  represent  the  front  7 (1)  as  the  level  set  (V>  = 0)  of  a function  <6.  To  motivate 
this  approach,  we  consider  the  example  of  an  expanding  circle.  Suppose  the  initial 
front  7 at  t = 0 is  a circle  in  the  xy-planc  (figure  3.1(a)).  We  imagine  that  the  circle 
is  the  level  set  {((.  = 0)  of  an  initial  surface  * = >fr(x,y,l  = 0)  in  /f3  (see  figure 
3.1(b)).  We  can  then  match  the  one-parameter  family  of  moving  curves  7(f)  with  a 
one-parameter  family  of  moving  surfaces  in  such  a way  that  the  level  set  (0  = 0) 
always  yields  the  moving  front  (figures  3.1(c)  & 3.1(d)). 

In  the  general  case,  let  7(0)  be  a closed,  noninlersecting,  (W  - 1)  dimensional 
hypersurface.  Let  V>(x,(),  x 6 RN , be  the  scalar  function  such  that  0(x,O)  = ±d(x), 
where  d(x)  is  the  signed  distance  from  x to  the  hypcrsurfacc  7(0).  We  use  the  plus 


Figure  3.1.  Level  set  formulation  of  equations  of  motion  - (a)  & (b)  show  the  curve  7 
and  the  surface  at  l = 0,  and  (c)  & (d)  show  the  curve  7 and  the  corresponding 
surface  t&(x,y)  at  timet. 

sign  if  x is  outside  7(0)  and  minus  sign  if  x is  inside.  Each  level  set  of  i/>  flows  along 
its  gradient  field  with  speed  f '(/().  The  gradient  V(f(x,f)  is  normal  to  the  (AT  - 1) 
dimensional  level  set  passing  through  x.  Now,  we  derive  the  equation  of  motion  for 
function  ip. 

Consider  the  motion  of  some  level  set  {i>  = 6’}.  Following  the  derivation  in  [30], 
let  x(l)  be  the  trajectory  of  a particle  located  on  this  level  set,  so 


The  particle  speed  dx/dt  in  the  direction  n normal  to  7(f)  is  given  by  the  speed 
function  F.  Thus, 

§rn=F-  (3-3) 

where  the  normal  vector  n is  given  by  n = Vi/'/  | Vf/|  |.  By  the  chain  rule  we  get, 


i(x(l),t)  = C. 


(3.2) 


V*  = 0 


(3.4) 


and  substitution  yields 


<h  + F | V0 1=  0, 


(3.5) 


with  an  initial  condition  0(x,O)  = ±d(x).  We  refer  to  equation  (3.5)  as  the  level 
set  “Hamilton-Jacobi"  formulation.  Note  that  at  any  time,  the  moving  front  -y(t)  is 
simply  the  level  set  (0(x,t)  = 0).  There  arc  several  advantages  to  this  approach. 
First,  since  the  underlying  coordinate  system  is  fixed,  discrete  mesh  points  used  in  the 
numerical  update  equations  do  not  move,  resulting  in  a stable  computation.  Topo- 
logical  changes  m the  front  can  be  handled  naturally  by  exploiting  the  property  that 
the  level  surface  (0  = 0)  need  not  be  simply  connected.  The  signed  distance  function 
0(x,i)  always  remains  a function,  even  if  the  level  surface  (0  = 0)  corresponding  to 
the  front  y(i)  changes  topology,  or  forms  sharp  corners.  The  geometric  and  differen- 
tial properties  of  7(l)  arc  captured  in  the  function  0 and  can  be  readily  extracted. 

As  an  example,  if  x € if*,  the  curvature  is  given  by 

K _ (iM1]- 20,0,0,,  + 0«0’) 

(02  + 0jj5/3  ■ (3.6) 

This  approach  can  also  be  easily  extended  to  higher  dimensions  and  appropriate 
expressions  can  be  obtained  for  the  mean  curvature  and  the  Gaussian  curvature  [30]. 
For  example,  the  mean  curvature  at  a point  (*,»,*)  on  a moving  surface  can  be 
computed  using  the  following  formula: 

A-  = + + *1)  + *M#  + 02)  - 20„0,0,  - 20„0,0,  - 20, ,0,0, 

(02  + 0J  + 02)3,!  - 

By  substituting  F(/( ) = 1 - eK  as  a typical  speed  function  in  equation  (3.5),  the 
equation  of  motion  becomes 


0i+  | V0  \=tK  | V0  | , 


Equation  (3.7)  resembles  a Hamillon-Jacobi  equation  with  viscosity,  where  “viscos- 
ity" refers  to  the  second-order  parabolic  right-hand  side.  This  equation  can  be  solved 
using  the  stable,  entropy-satisfying  finite  difference  schemes,  borrowed  from  the  lit- 
erature on  hyperbolic  conservation  laws  (see  [30]). 


SHAPE  RECOVERY  WITH  FRONT  PROPAGATION 


In  this  chapter,  we  describe  how  the  level  set  formulation  for  the  front  propagation 
problem  discussed  in  the  previous  chapter  can  be  used  for  shape  recovery.  There 


recovery.  In  the  former,  the  front  represents  a solid/liquid  interface  (crystal  growth) 
or  a boundary  separating  burnt  and  unburnt  regions  (dame  propagation).  In  these 
cases  the  computation  is  alive  as  long  as  there  remains  a physical  domain  into  which 
the  front  can  be  moved.  For  example,  the  flame  front  can  be  moved  as  long  as  there 
is  a region  to  be  burnt  and  it  hasn’t  crossed  the  physical  domain  in  which  the  solution 
is  sought.  On  the  contrary,  in  shape  recovery  the  front  represents  the  boundary  of 
an  evolving  shape.  Since  the  idea  is  to  extract  object  shapes  from  a given  image,  the 
front  should  be  forced  to  stop  in  the  vicinity  of  the  desired  object  boundaries.  This 
is  analogous  to  the  force  criterion  used  to  push  the  active  contour  model  towards 
desired  shapes.  We  define  the  final  shape  to  he  the  configuration  when  all  the  points 
on  the  front  come  to  a stop,  thereby  bringing  the  compulation  to  an  end. 

4J — SlgMing.Crilerion;  An  Image-based  Speed  Term 
Our  goal  now  is  to  define  a speed  function  from  the  image  data  that  can  be  applied 
on  the  propagating  front  as  a stopping  criterion.  In  general  the  function  F can  be 
split  into  two  components:  F = FA  + Fa.  The  term  Ft,  referred  to  as  the  advection 
term,  is  independent  of  the  moving  front’s  geometry.  The  front  uniformly  expands 
or  contracts  with  speed  FA  depending  on  its  sign  and  is  analogous  to  the  inflation 


force  defined  in  [9],  The  second  term  Fa,  is  the  part  which  depends  on  the  geometry 
of  the  front,  such  as  its  local  curvature.  This  (diffusion)  term  smooths  out  the  high 
curvature  regions  of  the  front  and  has  the  same  tegularizing  effect  on  the  front  as  the 
internal  deformation  energy  term  in  thin-platc-mcinbrane  splines  [17]  (see  the  figure 
(6.3)).  We  rewrite  equation  (3.7)  by  splitting  the  influence  of  F as 

* + ft  I V*  | +Fa  | Vd>  |=  0.  (4.1) 

First  consider  the  case  when  the  front  moves  with  a constant  speed,  i.c.  F = FA. 
To  this  if  we  add  a negative  speed  term  synthesised  front  the  image,  such  that  their 
sum  tends  to  aero  near  large  image  gradient  locations,  we  will  achieve  our  goal  of 
bringing  the  front  to  a slop  in  the  neighborhood  of  object  boundaries.  To  this  end, 
we  define  a negative  speed  Fi  to  be 

= W-*Uv)  11  VGV  * l(x'y)  I “"»)  • (4.2) 

where  M,  ard  hh  arc  the  maximum  and  minimum  values  of  the  magnitude  of  image 
gradient  | VG,  • l{x,y)  |,  (x,jj)  6 fi.  The  expression  G„  * / denotes  the  image 
convolved  with  a Gaussian  smoothing  filler  whose  characteristic  width  is  a.  Alter- 
nately, we  could  use  a smoothed  zero-crossing  image  to  synthesize  the  negative  speed 
function.  The  zero-crossing  image  is  produced  by  detecting  zero-crossings  in  the  func- 
tion v!6'„  * ;,  which  is  the  original  image  convolved  with  a Laplacian-of-Gaussian 
filter  whose  characteristic  width  is  <r.  The  equation  of  motion  with  the  addition  of 

*>  + (ft  + F,)  j Vi- 1=  0.  (4.3) 

F,  is  called  an  extension  of  F,  to  points  away  from  the  boundary  i((),  i.e.  at  points 
(*>»)  s (*>  - 7(f)).  and  is  equal  to  F,  on  i(l).  We  shall  return  to  the  issue  of 


extension  shortly.  The  value  of  F,  lies  in  the  range  [-F.,,0]  as  the  value  of  image 
gradient  varies  between  Afi  and  M,.  From  this  argument  it  is  clear  that  the  front 
gradually  attains  zero  speed  as  it  gets  closer  to  the  object  boundaries  and  eventually 

In  the  case  when  the  front  moves  with  a speed  that  is  a function  of  local  curvature, 
i.e.  Fq  yt  0,  it  is  not  possible  to  find  an  additive  speed  term  from  the  image  that  will 
cause  the  net  speed  of  the  front  to  approach  zero  in  the  neighborhood  of  a desired 
shape.  Instead,  we  multiply  the  speed  function  F = Fa  + Fa  with  a quantity 
The  term  k,,  which  is  defined  as 


has  values  that  are  closer  to  zero  in  regions  of  high  image  gradient  and  values  that 
are  closer  to  unity  in  regions  with  relatively  constant  intensity.  More  sophisticated 
stopping  criteria  can  be  synthesized  by  using  the  orientation  dependent  "steerable- 
filters  [16],  The  modified  equation  of  motion  is  given  by 


0<  + k,(Ft  + Fa)  | Vi/.  |=  0. 


We  now  come  to  an  important  juncture  in  our  discussion.  The  image-based  speed 
term,  be  it  F,  or  k,.  has  meaning  only  on  the  boundary  i(l),  i.c.  on  the  level  set 
{0  = 0}.  This  follows  from  the  fact  that  they  were  designed  to  force  the  propagating 
level  set  (0  = 0)  to  a complete  stop  in  the  neighborhood  of  an  object  boundary. 
However,  the  equation  of  motion  (4.3)  is  written  for  the  function  0,  which  is  made 
up  of  infinitely  many  level  curves.  In  other  words,  equations  (4.3)  k (4.5)  control 
the  evolution  of  a family  of  level  sets.  Therefore,  it  is  imperative  that  the  net  speed 
used  in  the  evolution  equation  has  a consistent  physical  meaning  for  all  the  level 


Figure  4.1.  Huygcn's  principle  construction 
sets,  i.e.  at  every  point  (*.*)  6 fl.  Speed  functions  such  as  Fa  which  are  functions 
of  geometric  properties  of  the  surface  = *(*,„),  can  be  readily  computed  at  any 
(*,»)  6 n.  However,  F,  is  not  such  a function.  It  derives  its  meaning  not  from  the 
geometry  of  0 but  from  the  configuration  of  the  level  set  {V  = 0}  in  the  image  plane. 
Thus,  our  goal  is  to  construct  an  image-based  speed  function  Fi  that  is  globally 
defined.  We  call  it  an  extension  of  F,  off  the  level  set  (0  = 0)  because  it  extends  the 
meaning  of  F,  to  other  level  sets.  Note  that  the  level  set  (0  = 0}  lies  in  the  image 
plane  and  therefore  F,  must  equal  F,  on  {0  = 0}.  The  same  argument  applies  to 
the  coefficient  b/. 

If  the  level  curves  arc  moving  with  a constant  speed,  i.e.  Fa  = 0,  then  at  any 
time  i,  a typical  level  set  (0  = C),  C 6 «,  is  a distance  C away  from  the  level  set 
{0  = 0)  (see  figure  4.1).  Observe  that  the  ahove  statement  is  a rephrased  version  of 
Huygen’s  principle  which,  from  a geometrical  standpoint,  stipulates  that  the  position 
of  a front  propagating  with  unit  speed  at  a given  lime  t should  consist  of  only  the 
set  or  points  located  a distance  t away  from  the  initial  front.  On  the  other  hand,  if 
Fa  # 0,  the  level  sets  will  violate  the  property  that  they  are  a constant  distance  away 
from  each  other.  However,  they  will  never  collide  and  cross  each  other  if  the  speed 


Figure  4.2.  Extension  of  image-based  speed  terms  to  other  level  sets 

function  F = Ft + Fa  is  continuous  (see  [14]).  With  the  above  remarks  in  mind  we 
slate  the  following: 

Fautrly  4- 9,1  Any  external  (image-basal)  speed  {unction  that  is  used  in  the  equation 
of  motion  written  for  the  function  0 should  not  cause  the  level  sets  to  collide  and  cross 
each  other  during  the  evolutionary  process. 

Recall  that  the  function  0(x,t)  has  been  initialised  to  d(x),  where  d(x)  is  the 
signed  distance  from  a point  X <=  SI  to  the  boundary  7(0).  Since  we  cannot  attribute 
any  geometric  meaning  to  the  function  F,  (*,)  at  points  away  from  the  level  set 
10  = 0),  we  look  for  a meaning  that  is  consistent  with  property  (4.2.1).  Therefore, 
the  question  to  ask  is:  what  is  the  value  of  F,  (or  hi)  at  a point  (x,  y)  lying  on  a level 
set  {0  = C)‘!  We  answer  this  question  in  the  following  construction  (sec  figure  4.2). 
Qlmlmcti"’)  4,gJ  The  value  of  F,  (in)  at  a point  P lying  on  a level  set  {0  = C) 
is  exactly  the  value  of  F,  (k,)  at  a point  Q,  such  that  point  Q is  a distance  C away 
from  P and  lies  on  the  level  set  (0=0). 

It  IS  easy  to  see  that  F,  reduces  to  F,  on  {0  = 0).  We  use  the  same  construction 
to  determine  the  value  of  hi  at  a point  (x,y)  lying  on  some  level  set  {0  = C). 


Note  that  if  th«  definition  of  a speed  function  adheres  to  construction  4.2.1,  then  it 
will  also  be  consistent  with  the  property  -1.2.1.  Thus,  having  ascribed  the  intent  of 
pseudodiffercntial  equations  (4.3)  il  (4.5)  in  the  context  of  shape  modeling,  wo  can 
use  finite  difference  schemes  to  solve  them  numerically.  Since  0 can  develop  comers 
and  sharp  gradients,  numerical  schemes  borrowed  from  hyperbolic  conservation  laws 
are  used  to  produce  stable  upwind  schemes.  Moreover,  the  equations  of  motion  can 
be  solved  on  a uniform  mesh  and  the  level  sets  can  be  moved  without  their  explicit 
construction. 

4.3  Shane  Recovery  in  .’ID 

In  this  section,  we  show  that  the  surface  shapes  of  three-dimensional  objects 
can  be  recovered  using  a scheme  that  is  analogous  to  the  one  described  above.  We 


level  set  {0  = 0}  of  a function  0(x,f),  x6ff.  The  function  0(x,f  = 0)  = ±d(x), 
where  d(x)  is  the  signed  distance  function  from  the  point  x to  the  initial  surface;  plus 
sign  is  used  if  x is  outside  the  initial  surface  and  minus  sign  if  x is  inside.  Each  level 
set  of  0 moves  along  its  gradient  field  with  speed  F(K).  Following  the  derivation  in 
chapter  2,  we  arrive  at  the  equation  of  motion  for  the  moving  surface,  namely 


0,  + F(K ) | V0  |=  0, 


(4.6) 


e.  Recall  that  the  surface 


at  any  time  I is  simply  the  level  set  (0  = 0). 

In  this  work,  we  only  address  the  issue  of  recovering  surface  shapes  from  3D 
images.  The  input  data  constitutes  a set  of  MRI  (Magnetic  Resonance  Imaging)  or 
CT  (Computed  Tomography)  images.  A stack  of  such  images  represents  an  image 
function  l(z,y,z).  Our  next  step  is  to  define  an  irnage-based  speed  term  which  will 


cause  llie  level  set  { l'  = 0}  to  slop  near  tlie  30  object  boundary.  To  this  end,  we 
apply  on  the  front  the  speed  equal  lo  the  product  of  f\l<)  and  a term  i/.  The  term 

has  values  closer  to  xero  in  the  vicinity  of  object  boundaries  and  values  that  arc  closer 
to  unity  in  regions  of  relatively  constant  image  intensity.  To  attribute  a consistent 

tension  to  ki  in  a manner  similar  to  the  construction  (4.2.1).  Therefore,  the  modified 
equation  of  motion  that  we  solve  to  extract  complex  surface  shapes  from  3D  images 


<h  + hHK)  I Vi>  |=  0. 


(4.8) 


NUMERICAL  SOLUTION 


The  equation  (3.7)  poses  an  initial  valued  problem.  It  is  rewritten  hero  as 


(jUi)  rn 


with  ’p{x,y,t  = 0)  = ± distance  from  (x,y)  to  q(0).  As  shown  in  Sethian  [38],  for 
£ > 0,  the  parabolic  right-hand  side  diffuses  sharp  gradients  and  forces  V*  to  stay 
smooth  at  all  values  of  1.  For  e = 0,  the  boundary  moves  with  unit  speed,  and  a 
comer  must  develop  from  smooth  initial  data.  Once  a corner  develops,  it  is  not  clear 


at  the  comer.  A variety  of  “weak”  solutions  which  propagate  the  curve  beyond  the 
occurrence  of  a singularity  are  possible.  Of  all  such  weak  solutions,  one  is  interested  in 
the  one  that  is  the  limit  of  smooth  solutions  as  e 0.  This  particular  weak  solution 
can  be  selected  with  the  help  of  a so-called  “entropy  condition",  see  |38|:  If  the  front 
is  viewed  as  a burning  Home,  then  once  a particle  is  burnt  it  slays  burnt.  Thus, 
approximations  to  the  spatial  derivative  are  sought  that  do  not  artificially  smooth 
sharp  cornets  and  which  pick  out  the  correct  entropy  solution  when  singularities 
develop.  The  schemes  given  in  (30,  39]  are  motivated  by  the  fact  that  the  entropy 
condition  for  propagating  fronts  is  identical  to  the  one  for  hyperbolic  conservation 
laws,  where  stable,  consistent,  entropy-satisfying  algorithms  have  a rich  history. 


In  this  section,  almost  without  a change,  we  present  the  arguments  discussed  in 
Sethian  [39],  For  complete  details  of  the  following  scheme,  we  refer  the  reader  to 
Osher  and  Sethian  [30,  39). 

First,  consider  the  ono-dimensional  version  of  the  level  set  equation,  with  constant 
normal  velocity  PA  = J.  Then  the  equation  (3.7)  becomes  a standard  liamilton- 
Jacobi  equation 

A + »(A)  = 0,  (5.2) 

where  U{4,,)  = -($) 1/2,  and  with  a given  initial  value  of  A Let  u = 4.  Taking  the 
derivative  with  respect  to  x,  equation  (S.2)  becomes 


u, + [«(«)[,  = 0,  (5.3) 

where  //(tr)  = — (u1)1^-.  Equation  (5.3)  is  a scalar  hyperbolic  conservation  law  in  one 
space  variable.  Solutions  can  develop  discontinuous  jumps,  known  as  shocks  even 

an  integral  version  of  the  conservation  law  which  admits  discontinuous  solutions  is 
studied.  Both  sides  of  equation  (5.3)  are  integrated  in  an  arbitrary  interval  [a,  6)  to 


sjf  = *[»(«,«)]  - ff[u(M)].  (5.4) 

u is  known  as  a tsea k solution  of  the  conservation  law  if  it  satisfies  the  above  integral 
equation.  Note  that  u need  not  be  differentiable  to  satisfy  the  integral  form  of  the 

When  will  a numerical  algorithm  approximate  the  correct,  entropy-satisfying  so- 
lution to  equation  (5.4)7  The  answer  is  found  in  the  following  definition: 


Vcp nilitn  9.1.1  ut  u,"  be  the  value  oj  u at  a mesh  paint  i Ax  at  lime  nAt.  / 
point  difference  scheme  is  said  to  be  in  conservation  form  if  there  exists  a f 
«(“ i.“j)  *«*  lllal  ‘he  scheme  can  be  written  in  the  form 
“"tl  ~ <&, ) ~ aW_, , «;') 


where  g(u,u)  = II  (u). 

This  definition  is  natural;  the  scheme  must  approximate  the  hyperbolic  conser- 
vation law,  subject  to  the  consistency  requirement  3(11,  u)  = II (u).  In  order  to 
guarantee  that  the  scheme  picks  out  the  correct  entropy-satisfying  weak  solution, 

ui-i»  11  < ’ ,UK^  u?+i-  fke  main  fact  is:  A conservative , monotone  scheme  produces  a 
solution  that  satisfies  the  entropy  condition.  Equation  (5.5)  is  a scheme  for  the  slope 
u,  which  must  be  converted  into  a scheme  for  i itself.  First  write  equation  (5.2)  with 
a forward  difference  in  lime  as, 


-AlH(u). 


(5.6) 


Since  the  numerical  flux  function  g approximates  II,  the  solution  to  equation  (5.6) 
may  be  approximated  by 

= ‘K  ~ A<s(“;-i/s,Ui+i/i) 

= A"  - AfjfOj’&.D+d,),  (5.7) 

where  g is  an  appropriate  numerical  flux  function  and  the  standard  definitions  of  the 
forward  and  the  backward  difference  operators  have  been  used,  namely, 


(5.8) 


Finally,  An  appropriate  numerical  flux  function  g is  required.  In  tire  special  case 
where  H(u)  may  be  written  as  a function  ofu2,  i.e.,  tf(u)  = /(u2)  for  some  function 
/,  one  can  use  the  Hamilton-Jacobi  flux  function  given  in  [30]: 

S(«.-I„,u, >./,)  = 

= /((roax(0;^,O))2  + (min(0;6O))2).  (5.9) 

This  conservative  monotone  scheme  is  upwind  method  in  that  it  differences  in  the 
direction  of  propagating  characteristics.  This  is  important,  since  it  imposes  bound- 
ary conditions  on  the  walls  of  a finite-sized  computational  box.  An  upwind  scheme 
automatically  differences  in  the  outward-flowing  direction  at  the  box  walls  if  the 
boundary  is  expanding,  thus  information  flows  out.  In  the  case  when  F*  = I,  so  that 
/(ua)  = — (u2)1/2,  equation  (5.7)  reduces  to 


*"*'  = *7  - At((max(O;d,0))J  + (min(0J*,O))2}''2.  (5.10) 
This  algorithm  produces  the  correct  entropy-satisfying  weak  solution  to  the  propa- 
gating front  problem. 

The  above  discussion  can  be  extended  to  problems  in  more  than  one  space  di- 
mension (sec  Osher  and  Scthiao  [30]).  Recall  that  the  original  intent  was  to  solve 
equations  (4.3)  and  (4.5).  In  two  dimensions,  the  scheme  given  in  equation  (5.10)  is 
extended  by  differencing  in  each  direction  to  produce  the  following  numerical  scheme 
for  equation  (4.3): 


Wj  - £l[Fx  + (Fr),  j]((max(O;d’iJ,0))a  + (min(O+*j,0))J 

+(max(Z)“i/',j,0))a  + (min(D+^1J-,0))J},/2. 


~ Cj  - A^sfi'jjMKmaxIOJiiijj.O))1  + (min(Oji/.,J,0))’  (5.12) 

+(max(O-0ij-,O))s  + (minfO+^.O))2)''2  - Al Fak,  | Vi- 1 . 


1 he  second  term  Fokt  I V0  | is  not  approximated  in  the  above  equatia 
use  a straightforward  central  difference  approximation  to  this  term. 


In  this  section  we  first  consider  a straightforward  approach  to  solve  the  equations 
(5.11)  & (5.13)  and  subsequently  show  that  the  time  complexity  of  our  algorithm  can 
be  improved  by  applying  the  update  procedure  to  a small  set  of  points  in  the  domain. 
In  other  words,  we  propagate  the  front  by  solving  the  governing  equation  at  a small 
number  of  points.  These  points  lie  in  a narrow  band  around  the  level  set  (0  = 0). 

It  should  be  observed  that  by  updating  the  level  set  function  on  a grid,  we  arc 
moving  the  level  sets  without  constructing  them  explicitly.  A particular  level  set 
(0  = 0),  which  models  the  shape  of  interest,  is  forced  to  a slop  in  the  vicinity  of 
desired  object  boundaries.  This  task  is  accomplished  by  applying  an  image-based 
speed  term  on  the  zero  set  (equation  (4.5)).  An  extension  function  is  constructed  to 
attribute  a consistent  meaning  to  the  image-based  speed  term  at  points  away  from 
the  zero  set  (see  discussion  in  section  3).  Therefore,  a straightforward  algorithm  con- 
sists of  advancing  from  one  time  step  to  the  next  as  follows: 

Algorithm  1 

1.  At  each  grid  point  (lAx.jAy),  where  Ax  and  Ay  are  step  sizes  in  cither  coor- 
dinate directions,  the  extension  of  image-based  speed  term  is  computed.  This 
is  done  in  accordance  with  the  construction  described  in  previous  section;  ix., 
by  searching  for  a point  q which  lies  on  the  lovol  set  (0  = 0),  and  is  closest 
to  the  point  (lAx.jAy).  The  value  of  image-based  speed  term  at  the  current 
point  is  simply  its  value  at  the  point  q. 


2.  With  the  value  of  extended  speed  term  (£y),v,  and  calculate  using 
the  upwind,  finite  difference  scheme  given  in  equation  (5.13). 

3.  Construct  an  approximation  for  the  level  set  {(f  = 0)  from  $"H.  This  is 
required  to  visualise  the  current  position  of  the  front  in  the  image  plane.  A 
piecewise  linear  approximation  for  the  front  y(l)  is  constructed  as  follows.  Given 
a cell  C(i,j),  ifmaxW,j,*(1j,*J>1,v.1+1Jtl)  < 0 or  min(*J,0i+1J,AJ+lf 
&+ij+i)  > 0,  then  C(i,j)  f 7(1)  and  is  ignored,  else,  the  entry  and  exit  points 
where  V = 0 are  found  by  linear  interpolation.  This  provides  two  nodes  on  7(<) 
and  thus,  one  of  the  line  segments  which  form  the  approximation  to  7(t),  The 
collection  of  all  such  line  segments  constitutes  the  approximation  to  the  level 
set  (0  = 0),  which  is  used  for  future  evaluation  of  the  image-based  speed  term 
in  the  update  equation  (5.13). 

4.  Replace  n by  n 4- 1 and  return  to  step  1 , 

It  is  easy  to  see  from  the  above  algorithm  that  the  most  expensive  step  is  the 
computation  of  the  extension  for  image-based  speed  term.  This  is  because  at  each  grid 
point,  we  must  search  for  the  closest  point  lying  on  the  level  set  (0  = 0).  Moreover, 
if  Fa  = 0,  then  the  stability  requirement  for  the  explicit  method  for  solving  our  level 
set  equation  is  AI  = O(Ax).  For  the  full  equation  (5.13),  the  stability  requirement 
is  At  = O(Ax’).  This  could  potentially  force  a very  small  time  step  for  fine  grids. 
These  two  effects,  individually  and  compounded,  make  the  computation  exceedingly 
slow.  A marginal  improvement  in  performance  can  be  achieved  by  evaluating  the 
extension  only  once  every  k iterations.  In  other  words,  we  take  k steps  in  time  without 
recomputing  the  image-based  speed  (*,)y.  Alternately,  we  could  down-sample  the 
image  and  perform  our  calculations  at  a lower  resolution.  On  the  other  hand,  we 
run  the  risk  of  losing  accuracy  while  computing  on  a coarse  grid.  One  approach  is  to 


Figure  5.1.  A narrow-band  of  width  5 around  tiie  level  set  { u = OJ. 
embed  our  level  set  algorithm  in  a multiresolution  framework.  Instead,  we  improve 
the  time  complexity  of  our  algorithm  by  updating  the  level  set  function  y only  at  a 
small  set  of  points  on  the  grid. 

The  basic  idea  behind  this  new  scheme  stems  from  the  observation  that  the  front 
can  be  moved  by  updating  the  level  set  function  at  a small  set  of  points  in  the 
neighborhood  of  the  aero  set  instead  of  updating  it  at  all  the  points  on  the  grid.  In 
figure  (5.1)  the  bold  curve  depicts  the  level  set  {0  = 0)  and  the  shaded  region  around 
it  is  the  narrow-band.  The  narrow-band  is  bounded  on  cither  side  by  two  curves  which 
are  a distance  4 apart,  i.c.,  the  two  curves  are  the  level  sets  {(4  = ±4/2}.  The  value 
of  4 determines  the  number  of  grid  points  that  fall  within  the  narrow-band.  Since, 
during  a given  time  step  the  value  of  0i;  is  not  updated  at  points  lying  outside  the 
narrow-band,  the  level  sets  (0  >|  4/2  |}  remain  stationary.  The  zero  set  which  lies 
inside  moves  until  it  collides  with  the  boundary  of  the  narrow-band.  Which  boundary 
the  front  collides  with  depends  on  whether  it  is  propagating  inward  or  outward;  either 
which  way,  it  cannot  move  past  the  narrow-band's  boundary. 


Therefore,  os  o consequence  of  our  update  strategy,  the  front  can  be  moved 
through  a maximum  distance  of  d/2,  cither  inward  or  outward,  at  which  point  we 
must  rebuild  an  appropriate  narrow  band.  We  reinitialize  the  t/r  function  by  treating 
the  current  zero  sot  configuration,  i.e.,  (it  = 0},  as  the  initial  curve  7(0)  (see  section 
2),  Note  that  the  reinitialization  procedure  must  account  for  the  case  when  (V1  = 0} 
changes  topology.  This  procedure  will  restore  the  meaning  of  ^ function  by  correcting 
the  inaccuracies  introduced  as  a result  of  our  update  algorithm.  Once  a new  i/>  func- 
tion is  defined  on  the  grid,  we  can  create  a new  narrow-band  around  the  zero  set,  and 
go  through  another  set  of,  say  /,  iterations  in  time  to  move  the  front  ahead  by  a dis- 
tance equal  to  d/2.  The  value  of  / is  set  to  the  number  of  time  steps  required  to  move 
the  front  by  a distance  roughly  equal  to  d/2.  This  choice  depends  on  some  experi- 
mentation. Thus,  a faster  algorithm  for  shape  recovery  consists  of  the  following  steps: 


Algorithm  2 

1.  Set  the  iteration  number  m = 0 and  go  to  step  2. 

2.  At  each  grid  point  (i,  j)  lying  inside  the  narrow-band,  compute  the  extension 
hi  of  image-based  speed  term. 

3.  With  the  above  value  of  extended  speed  term  (if)y,  and  calculate  $".+* 
using  the  upwind,  finite  difference  scheme  given  in  equation  (5.13). 

4.  Construct  a polygonal  approximation  for  the  level  set  = 0}  from  A 
contour  tracing  procedure  is  used  to  obtain  a polygonal  approximation.  Given 
a cell  (i,j)  which  contains  7(f),  this  procedure  traces  the  contour  by  scanning 
the  neighboring  cells  in  order  to  find  the  next  cell  which  contains  7(1).  Once 
such  a cell  is  found,  the  process  is  repeated  until  the  contour  closes  on  itself. 


The  set  of  nodes  visiled  during  this  tracing  process  constitutes  the  polygonal 
approximation  to  7(t).  In  general,  to  collect  all  the  closed  contours,  the  above 
tracing  procedure  is  started  at  a new,  as  yet  unvisitod,  cell  which  contains 
the  level  set  {V  = 0}.  A polygonal  approximation  is  required  in  step  2 for 
the  evaluation  of  image-based  speed  term  and  more  importantly,  in  step  6 for 
reinitializing  the  ij  function. 

5.  Increment  m by  one.  If  the  value  of  m equals  I,  go  to  step  6,  else,  go  to  step  2. 

6.  Compute  the  value  of  signed  distance  function  * by  treating  the  polygonal 
approximation  of  {0  = 0}  as  the  initial  contour  7<0).  As  mentioned  earlier, 
a more  general  method  of  reinitialization  is  required  when  {gf>  = 0}  changes 
topology.  Go  to  step  1. 

The  issue  of  boundary  conditions  needs  to  be  addressed  at  this  time.  In  our 
original  method,  the  level  set  function  * was  initialized  and  updated  on  a square 
grid.  To  update  at  a boundary  point  (i.j),  it  is  necessary  to  specify  the  boundary 
conditions  in  order  to  numerically  approximate  the  derivatives.  We  have  chosen 
free-end  boundary  conditio™,  i.e„  if  a derivative  computation  at  a given  boundary 
point  (i,i)  involves  a point  that  falls  outside  the  computational  domain,  then  that 
derivative  is  set  to  zero.  In  our  new  approach,  since  we  only  update  * at  points 
lying  in  the  narrow-band,  the  issue  of  specifying  boundary  conditions  for  points 
lying  on  the  edge  of  the  band  becomes  pertinent.  The  focus  of  front  propagation 
in  the  context  of  shape  recovery  is  merely  rapid  advancement  of  the  front.  Therefore, 
at  the  expense  of  a little  inaccuracy  in  the  zero  set  configuration,  we  ignore  the 
sensitive  issue  of  specifying  correct  boundary  conditions.  In  our  implementation, 
derivatives  are  evaluated  normally  even  if  it  involves  a point  that  doe.  not  belong 


in  tile  narrow-band.  On  the  oilier  hand,  in  applications  such  as  crystal  growth,  and 
flame  propagation,  accurate  specification  of  boundary  conditions  is  imperative. 

We  now  show  that  this  new  faster  approach  provides  a correct  approximation  to 
the  propagating  front  problem.  In  figure  (5.2),  we  show  the  result  of  applying  narrow- 
band  updation  algorithm  to  a star  shaped  front  propagating  with  speed  F = -K, 
where  K is  the  curvature  as  in  equation  (3.6).  The  calculation  was  done  on  a unit 
box  with  64  points  in  either  direction,  and  a time  step  of  At  = 0.0003  was  employed. 
The  width  of  the  narrow  band  has  been  set  to  6 = 0.075,  and  the  V>  function  was 
recomputed  once  every  (f  =)  40  time  steps.  In  figure  5.2(a),  we  show  the  initial 
curve  along  with  the  level  sets  (i/i  = ±i0.04),  where  i e [1..5].  After  40  narrow-band 
updations  (figure  5.2(b)),  only  the  level  sets  {<!>  <|  0.08  |)  move  and  the  rest  remain 
stationary.  We  note  the  inconsistency  between  the  level  sets  lying  on  either  side  of 
the  narrow-band,  making  the  reinitialisation  step  necessary  in  order  to  restore  the 
meaning  of  the  i/>  function.  Following  the  reinitialisation  step,  another  40  update 
steps  are  applied  (figure  5.2(c)),  which  “diffuses”  the  high  curvature  regions  of  the 
front  even  further.  In  subsequent  figures,  the  results  of  repeatedly  applying  the 
same  strategy  are  shown.  Finally,  in  figure  5.2(f),  the  peaks  and  troughs  on  the 
front  get  completely  diffused,  and  it  attains  a smooth  circular  configuration  after  4 
reinitialisation  steps  and  a total  of  200  time  steps. 

We  conclude  this  section  by  highlighting  some  advantages  of  the  narrow-band 
updation  scheme.  The  first  is  the  obvious  improvement  in  time  complexity  as  a 
result  of  considering  only  a fraction  or  the  total  number  of  grid  points  for  updation. 
The  second  is  related  to  the  issue  of  computing  the  value  of  image-based  speed  term 
at  points  away  from  the  scro  set,  i.e.,  the  extension.  In  our  original  method,  the 
extension  of  image-based  speed  term  could  be  unreliable  at  points  that  are  relatively 
far  from  the  zero  sot.  In  addition,  when  the  front  is  a circular  one,  it  is  not  clear 


(e)  After  80  iterations 


(d)  After  120  iterations 


(e)  Alter  160  iterations  (f)  Alter  200  iterations 


Figure  5.2.  Narrow-band  updation  algorithm  applied  to  a star-shaped  front 
gating  with  speed  F = —K.  Calculations  were  done  on  a «1  x (il  grid  with 
step  Al  = 0.0003.  V’  was  recomputed  after  every  '10  time  iterations. 


»l.at  value  of  image-based  speed  term  to  use  at  its  center.  Here,  the  fact  that  the 
center  of  a circle  is  equidistant  from  all  the  points  lying  on  it  causes  an  ambiguity  in 
the  choice  of  the  shortest  distance  point  on  the  circular  front.  Clearly,  this  effect  is 
undesirable.  The  narrow-band  updation  algorithm  does  not  require  the  image-based 
speed  term's  evaluation  at  far  away  points.  Therefore,  the  extension  is  very  accurate 
and  reliable. 

5,3  Numerical  Issues  in  SI) 

In  this  section,  we  provide  some  details  of  the  numerical  method  based  upon  the 
Eulerian  formulation  [30)  that  is  used  to  solve  the  equation  (4.8).  To  begin,  the 
discussion  in  the  previous  sections  can  be  readily  extended  to  3D  by  differencing  in 
each  space  direction  to  produce  the  following  numerical  scheme  to  approximate  the 
surface  motion: 

0))J  + (mln(DJ  0))J  (5.13) 

+(max(D-^,t,0))J  + (min(OJ>i,M,(l))J 
+(max(D;*J>,0))’  + (min(D»(|i1J,J,0)),|,/1  - At(eK)(k,)  | VV  I . 
Note  that  the  level  surfaces  are  moving  with  speed  F = i,(l  - eK),  where  K is  the 

In  the  context  of  of  3D  shape  recovery,  we  move  the  level  surfaces  by  updating  the 
function  on  a three  dimensional  grid  according  to  the  above  update  equation.  The 
level  surface  {0  = 0),  which  models  the  shape  of  interest,  is  forced  to  a stop  in  the 
vicinity  of  desired  object  boundaries  by  applying  an  image- based  speed  coefficient  on 
it.  In  order  to  to  compute  the  value  of  image-based  speed  coefficient  at  points  away 
from  the  zero  set,  an  extension  is  constructed.  More  specifically,  the  shape  recovery 
procedure  in  3D  consists  of  the  following  stops: 


Algorithm  3 


1.  Al  each  grid  point  (t'Ax,jA»,AAs),  the  extension  of  image-based  speed  coclli- 

2.  With  the  value  of  extended  speed  coefficient  and  Wjjn  calculate 

using  the  upwind,  finite  difference  scheme  given  in  oquation  (5.14). 

3.  Construct  an  approximation  for  the  level  surface  (0  = 0}  from  0"^.  This  is 
done  in  the  following  fashion.  Across  every  two  vertices  of  a given  cell  C(i,j,  fc), 
the  value  of  'll  function  is  sampled.  If  the  0 function  changes  sign  across  a given 
edge,  the  exact  point  at  which  (S  = 0 is  found  by  linear  interpolation.  Note  that 
at  each  cell,  the  0 function  can  change  sign  across  a number  of  edges  that  lies 


to  the  level  surface  {</>  = 0),  which  is  used  for  the  evaluation  of  k,  during  the 
next  time  step. 

4.  Replace  rt  by  n + 1 and  return  to  step  1. 

In  the  previous  section,  it  has  been  established  that  the  extension  calculation 
is  very  expensive.  This  observation  has  motivated  the  narrow-band  approach  for 
front  propagation.  The  situation  in  3D  is  even  worse  due  to  an  extra  dimension.  In 
addition,  there  is  a large  degree  of  repetition  in  the  sero  set  due  to  the  lack  of  an 
organized  method  in  which  to  populate  the  zero  set.  By  organized  method  we  mean  a 
3D  analog  or  a 2D  contour  tracing  procedure.  This  limitation  adds  to  the  size  of  the 
zero  set  and  subsequently  to  the  time  complexity  of  extension  calculation.  In  2D,  we 
alleviated  this  situation  by  adopting  a narrow-band  update  strategy,  but  implement- 
ing a similar  scheme  in  3D  could  prove  overly  cumbersome.  Therefore,  we  exploit 


the  inherent  data  independence  that  is  present  in  our  grid  calculation,  and  devise  a 
parallel  implementation  to  move  the  level  surfaces.  The  parallel  implementation  was 
done  on  a KSR  (Kendall  Square  Research)  shared  memory  parallel  machine.  For  a 
grid  with  A'3  points,  we  have  used  N processors  and  realised  a substantial  improve- 
ment in  the  running  time  of  our  3D  shape  recovery  algorithm.  Steps  1 &-  2 of  the 
above  algorithm  are  done  in  parallel  and  step  3 is  executed  in  serial  inode. 


CHAPTER  6 

EXPERIMENTAL  RESULTS 

In  this  chapter  we  present  several  shape  recovery  results  that  were  obtained  by  ap- 
plying the  narrow-hand  level  set  algorithm  to  image  data.  Subsequently,  we  present 
some  results  in  3D.  First,  in  2D,  given  an  image,  our  method  requires  the  user  to 
provide  an  initial  contour  -y(0).  As  we  shall  sec,  there  is  absolutely  no  restriction  on 
where  the  initial  contour  can  be  placed  in  the  image  plane  as  long  as  it  is  inside  a 
desired  shape  or  encloses  all  the  constituent  shapes.  This  feature  is  of  paramount 
importance  in  the  context  of  automatic  shape  recovery.  Our  front  seeks  the  object 
boundaries  by  either  propagating  outward  in  the  normal  direction  or  propagating 
inward  in  the  negative  norma]  direction.  This  choice  is  made  at  the  lime  of  initial- 
ization. Note  that  after  the  specification  of  initial  shape  of  7(0),  our  algorithm  docs 
not  require  any  further  user  interaction. 

The  initial  value  of  the  function  i.c.,  <|i(x,0)  is  computed  from  7(0).  We  first 
discretize  the  level  set  function  tfi  on  the  image  plane  and  denote  <hj  as  the  value  of 
^ at  a grid  point  (iAz.jAy),  where  Az  and  Ay  are  step  sizes  in  either  coordinate 
directions.  In  our  implementation,  since  we  usually  work  with  2*  x 2*  images,  the 
computational  domain  is  a square  with  Az  = Ay  = A.  We  define  the  distance  from 
a point  (i,j)  to  the  initial  curvo  to  be  the  shortest  distance  from  (i,j)  to  7(0).  The 
magnitude  of  is  set  to  this  value.  We  use  the  plus  sign  if  (i,j)  is  outside  7(0) 
and  minus  sign  if  (1,  j)  is  inside.  Once  the  value  of  i/1,.,  is  computed  at  time  1 = 0 
by  following  the  above  procedure,  we  use  the  update  equations  from  the  previous 
chapter  to  move  the  front. 
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6.1  2D 


-Recovery  Results 


We  now  present  our  shape  recovery  results  in  2D.  We  first  consider  a 256  x 256  CT 
(computed  tomography)  image  of  an  abdominal  section  shown  in  figure  6.1(a).  Our 
intention  is  to  recover  the  shape  of  the  liver  in  this  particular  slice.  The  function  ip  has 
been  discretixed  on  a 128  x 128  mesh,  i.e.,  calculations  are  performed  at  every  second 
pixel.  In  figure  6.2(a),  we  show  the  closed  contour  that  the  user  places  inside  the 
desired  shape  at  lime  t = 0.  The  function  ip  is  then  made  to  propagate  in  the  normal 
direction  with  speed  F = h/(— 1.0  — 0.025/C).  We  have  employed  the  narrow-band 
updation  algorithm  to  move  the  front.  The  time  step  size  was  sot  to  At  = 0,0005 
and  the  ip  function  was  recomputed  after  every  50  time  steps.  Figure  6.1(b)  shows 
the  image-based  speed  term  which  is  synthesised  according  to  equation  (4.4).  This 
term  when  applied  on  the  propagating  front,  acts  as  a stopping  criterion,  thereby 
causing  the  front  to  come  to  a rest  in  the  neighborhood  of  a desired  shape.  Note 
that  in  figure  6.1(b),  fc/(z,!f)  values  lying  in  the  interval  [0..1]  have  been  mapped  into 
the  interval  [0..255].  In  figures  6.2(b)  through  6.2(c)  we  depict  the  configuration  of 
the  level  set  (C  — 0}  at  four  intermediate  time  instants.  The  final  result  is  achieved 
after  575  time  iterations  and  is  shown  in  figure  6.2(f).  It  should  be  noted  that  our 
method  does  not  require  that  the  initial  contour  be  placed  close  to  the  object  boundary 
In  addition,  observe  how  the  front  overshoots  all  the  isolated  spurious  edges  present 
inside  the  shape  (see  figure  6.1(b))  and  settles  in  the  neighborhood  of  edges  which 
correspond  to  the  true  shape.  This  feature  is  a consequence  of  eK  component  in  the 
speed  which  diffuses  regions  of  high  curvature  on  the  front  and  forces  it  to  attain  a 
smooth  shape. 

As  mentioned  in  section  3,  smoothness  of  the  front  can  be  controlled  by  choosing 
an  appropriate  curvature  component  in  the  speed  function  F = 1 -eK.  The  objective 


of  our  next  experiment  is  to  demonstrate  smoothness  control  in  the  context  of  shape 
recovery.  In  figures  6.3(a)  through  6.3(c),  we  show  the  results  of  applying  our  narrow- 
baud  shape  recovery  algorithm  to  an  image  consisting  of  three  synthetic  shapes. 
Initialisation  was  performed  by  drawing  a curve  enclosing  each  one  of  the  three 
shapes.  We  compute  the  signed  distance  function  gfr  (*,  y)  from  these  curves.  The  level 
sets  arc  then  made  to  propagate  with  speed  F = i, (1.0 -eK).  First,  as  shown  in  figure 
6.3(a),  we  perform  shape  recovery  with  the  value  of  e = 0.05.  The  process  is  repeated 
with  different  values  of  e;  0.25  in  figure  6.3(b)  and  0.75  in  figure  6.3(c).  Clearly,  with 
every  increment  in  the  value  of  e,  the  level  set  (f  = 0}  attains  a configuration  that 
is  relatively  smoother.  This  is  analogous  to  the  controlled-continuity  provided  by 
thin-plate-membrane  splines  (41,  17). 

In  our  third  experiment  we  recover  the  complicated  structure  of  an  arterial  tree. 
The  real  image  has  been  obtained  by  clipping  a portion  of  a digital  subtraction 
angiogram.  This  is  an  example  of  a shape  with  brandies  and  significant  protrusions. 
In  this  experiment  we  compare  the  performance  of  our  sdicinc  with  the  active  contour 


Figure  6.2.  Recovery  of  the  liver  shape  from  a CT  image  of  an  abdominal  section 
Narrow-band  compulation  was  done  on  a 128  x 128  grid  the  front  was  marie  to 
propagate  with  speed  F = *,(-1.0  - 0.025A')  and  the  time  step  At  was  set  to 
0.0005.  til  was  recomputed  once  every  50  time  steps. 


(a)  £ = 0.05  (h)  c = 0.25 


Figure  6-3-  Smoothness  control  in  shape  recovery  can  be  achieved  by  varying  the 
curvature  component  in  the  speed  /■'  = £,(1.0  - cK). 

model  to  bring  the  limitations  of  the  later  into  focus.  We  first  attempt  to  reconstruct 
the  complex  arterial  structure  using  a snake  model  with  inflation  forces  [9].  In  figures 

6.4(a)  through  6.4(i),  we  show  a sequence  of  pictures  depicting  the  snake  configuration 

in  the  image.  We  present  the  final  equilibrium  state  of  the  snake  in  figures  6.4(c), 
6.4(f),  & 6.4(i)  corresponding  to  three  distinct  initialisations,  one  better  than  the 


barely  recovers  the  main  stem  of  the  artery  and  completely  fails  to  account  for  the 
brandies.  Two  prominent  limitations  of  the  snake  model  immediately  come  into  light. 
The  first  is  the  dependence  of  final  result  on  the  initial  configuration.  This  is  due  to 
the  existence  of  multiple  local  minima  in  the  (nonconvex)  energy  functional  which 
the  numerical  procedure  explicitly  minimizes.  The  second  feature  is  the  inability  of 
snake  model  to  attain  a stable  shape  with  protrusions.  Observe  how  in  the  third  case, 
despite  a good  initialization  (figure  6.4(g)),  the  snake  “snaps"  back  into  a relatively 
“bumpless”  configuration  in  figure  6.4(h).  This  inadequacy  stems  from  snake's  arc- 
length  (elasticity)  and  curvature  (rigidity)  minimizing  nature.  Snake  prefers  regular 
shapes  because  shapes  with  protrusions  have  very  high  deformation  energies.  Now, 


wo  apply  our  level  set  algorithm  to  reconstruct  the  same  shape.  After  initialization  in 
figure  6.5(a),  the  front  is  made  to  propagate  iu  the  normal  direction.  We  employ  the 
narrow-band  updation  scheme  with  a band  width  of  6 = 0,075  to  move  tile  front.  It 
can  be  seen  that  in  subsequent  frames  the  front  literally  "flows"  into  the  brandies  and 
finally  in  6.5(f)  it  completely  reconstructs  the  complex  tree  structure.  The  advantages 
of  our  scheme  are  quite  apparent  from  this  example.  Since  our  front  advancement 
process  docs  not  involve  optimization  of  any  quantity,  the  shape  recovery  results  we 
obtain  are  independent  of  initialization.  In  addition,  a single  instance  of  our  shape 
model  "sprouts”  branches  and  recovers  all  the  connected  components  of  a given  shape. 
In  figures  6.6(a)-(i)  we  plot  the  other  level  sets  to  elucidate  the  process  of  extending 
the  image-based  speed  function  to  points  away  from  the  zero  set.  AH  calculations 
were  carried  out  on  a 64  x 64  grid  and  a time  step  At  = 0.001  was  used. 

In  the  next  experiment,  we  depict  a situation  when  the  front  undergoes  a topo- 
logical transformation  to  reconstruct  the  constituent  shapes  in  an  image.  The  image 
shown  in  figure  6.7(a)  consists  of  three  distinct  shapes.  Initial  curve  is  placed  in  such 
a way  that  it  envelops  all  the  objects.  The  front  is  then  advanced  in  the  direction  of 
negative  normal.  Alternately,  we  could  perform  the  initialization  by  placing  a curve 
in  each  one  of  the  individual  shapes  and  propagating  them  in  the  normal  direction. 
We  choose  the  former  option.  The  level  set  = 0)  first  wraps  itself  tightly  around 
the  objects  (see  figures  6.7(c)  & 6.7(d))  and  subsequently  splits  into  four  separate 
closed  curves  (figure  6.7(e)).  While  the  first  three  closed  segments  of  {if  = 0)  recover 
the  three  distinct  shapes,  the  one  in  the  middle  (see  figure  6.7(c)),  since  it  does  not 
enclose  any  object,  eventually  disappears.  Figure  6.7(f)  shows  the  final  result.  Again 
it  should  be  noted  that  a single  instance  of  our  shape  model  dynamically  splits  into 
three  instances  to  represent  each  object.  To  show  the  working  of  our  algorithm,  in 


figures  G.8(a)-(1)  wo  show  the  lovol  sets  (0  = ±iC),  i e (-5..  + 5],  with  C = 0.05. 
The  function  v was  discretized  on  a 64  x 64  grid  and  At  was  set  to  0.001. 

Lastly,  we  show  that  our  approach  can  also  be  used  to  recover  shapes  with  holes. 
We  do  this  by  applying  our  method  to  extract  the  shapes  of  hand-written  characters. 
The  characters  “A"  and  "B”  in  the  figure  (6.9)  arc  examples  of  shapes  with  holes. 
Recovery  of  hand-written  character  shapes  finds  application  in  the  area  of  character 
recognition.  These  shapes  can  also  serve  as  an  input  stage  to  algorithms  computing 
medial  axis  transform  [21,  8]  or  skeletons  [29,  23).  The  attractive  feature  here  is 

recovered  without  separate  initialisations.  In  figure  6.9(a),  we  show  the  initial  contour 
which  encloses  ail  the  characters.  This  contour  is  then  made  to  propagate  inward  with 
a constant  speed.  Figures  6.9(b)  shows  an  intermediate  stage  in  the  front's  evolution 
and  in  6.9(c),  the  front  splits  into  four  separate  contours.  The  calculation  comes  to  a 
hall  when  in  figure  6.9(d),  the  level  set  {0  = 0)  recovers  the  outer  boundaries  of  four 
separate  characters.  Unlike  the  characters  “A"  & *B"  for  which  we  need  to  extract 
the  inner  boundary  for  their  complete  shape  description,  for  the  characters  “S"  & 
“H",  the  outer  boundary  completely  describes  their  shape.  In  the  second  stage  of  our 
computation,  wo  treat  tile  zero  set  configuration  in  figure  6.9(d)  as  an  initial  state, 
and  propagate  all  the  fronts  inward  by  momentarily  relaxing  the  image-based  speed 
term.  This  causes  the  zero  set  to  move  iuto  the  character  shapes  as  shown  in  figure 
6.9(e).  Finally  in  figure  6.9(f),  the  zero  set  comes  to  stop  in  the  neighborhood  of  holes 
that  are  present  in  the  shapes  of  characters  “A"  U “B”,  thereby  achieving  a complete 
shape  recovery.  The  calculations  for  this  experiment  were  done  on  a 128  X 128  grid 
and  the  time  step  At  was  set  to  0.00025. 


In  this  section,  we  present  some  front  propagation  and  shape  recovery  examples 
in  3D.  The  first  two  experiments  demonstrate  the  basic  surface  motion  scheme  in  3D 
with  a constant  and  curvature  dependent  speed.  To  visualise  tile  level  set  (v'  = 0),  we 
use  a telrahcdronization  algorithm.  Given  a function  on  a grid,  this  algorithm 
first  marks  off  all  the  points  with  negative  V'  values  (here  we  assume  that  the  points 
lying  inside  the  surface  have  negative  V'  values  and  the  points  lying  outside  have 
positive  v'  values).  A Lcirahcdrouixaliou  procedure  is  then  applied  to  these  points. 
The  triangular  facets  that  lie  on  the  boundary  of  the  tetrahedron  constitutes  our 
approximation  to  the  level  surface  {</i  = 0).  Using  these  surface  triangles  the  level 
set  (0  = 0),  which  forms  our  3D  surface  model,  is  rendered  as  a shaded  shell. 

In  our  first  3D  experiment,  we  evolve  a toroidal  surface  with  constant  speed  F — 1 . 
The  computational  box  is  a cubical  one  which  extends  from  -0.75  to  0.75  in  all  three 
coordinate  directions.  Initial  value  of  (&  is  defined  according  to  the  following  equation: 


The  set  of  all  points  (z,jf,z)  at  which  the  value  of  gf>  = 0,  lie  on  the  surface  of  the 
torus  (see  figure  6.10(a)).  We  discretize  the  tp  on  a 32  x 32  x 32  grid  and  follow  a 
nested  set  of  toroidal  shapes  by  updating  it  on  the  grid  according  to  the  equation 
(5.14).  Figures  6.10(b)-6.10(d)  depict  three  intermediate  stages  wherein  the  inner 
radius  gradually  becomes  smaller  until  in  figure  6.10(e)  the  inner  radius  collapses 
to  zero,  normals  collide,  and  topology  changes.  The  surface  at  this  stage  looks  like 
a dented  sphere  and  will  asymptotically  attain  a spherical  configuration.  In  figure 
6.10(f),  the  level  set  (d*  = 0)  intersects  with  the  edge  of  the  computational  box 
causing  the  slicing  of  the  toroidal  surface. 


In  our  next  experiment,  we  collapse  a “pinched"  superquadric  shape  under  its 
mean  curvature.  The  initialisation  was  done  using  the  implicit  function  form  of  a 
superquadric,  namely 

- [{(ft)*  * &a*r  * <ar]  - iji  •• 

where  t|,tj  > 0 ore  “squareness”  parameters.  By  varying  the  values  of  <,  and  c2, 
several  interesting  shapes  can  be  generated.  Flat  beveled  shapes  are  produced  when 
either  £|  or  tj  = 2 and  pinched  shapes  are  produced  when  cither  <i  or  cj  > 2.  Figure 
6.11(a)  shows  the  initial  superquadric  shape  with  sharp  corners.  This  shape  is  made 
to  move  with  speed  F = —K,  where  K is  the  mean  curvature.  As  shown  in  figures 
6.11(b)-  6.11(d),  the  high  curvature  regions  corresponding  to  the  sharp  corners  on 
the  surface  smoothly  diffuse  until  the  superquadric  collapses  into  an  ellipsoid  in  figure 
6.11(c).  Note  that  an  ellipsoid  is  a superquadric  with  c,  = t2  = 1.  Since  the  mean 
curvature  on  the  surface  of  an  ellipsoid  is  constant,  the  surface  eventually  collapses 
to  a point  without  changing  its  shape  (see  the  smaller  ellipsoid  in  figure  6.11(f)). 
Calculations  were  done  on  a 32  x 32  x 32  grid  with  a time  stop  A t = 0.0001. 

In  our  last  experiment,  we  recover  the  shape  of  a flat  superquadric  using  the  level 
set  front  propagation  scheme  in  3D.  Image  data  for  this  experiment  consists  of  32 
images  each  with  a particular  cross  section  of  the  superquadric.  The  image-based 
speed  term  4/  is  computed  from  these  images  according  to  equation  (4.7).  A sphere, 
which  is  the  level  surface  {0  = 0}  of  a function  = x1  + y*  + t2  - 0.01, 
forms  our  initialization  (sec  figure  6.12(a)).  This  initial  surface  is  moved  with  speed 
F = i,  by  updating  the  value  of  ^ on  a discrete  grid  according  to  the  details  given  in 
algorithm  3.  The  initial  surface  expands  smoothly  in  all  directions  until  a portion  of  it 
collides  with  the  superquadric  boundary.  At  points  with  high  gradient,  the  k,  values 
are  close  to  aero  and  cause  the  zero  set  to  locally  come  to  stop  near  the  boundary 


n is  depicted  in  figures  6.12(b)-  6.12(e), 
wherein  the  initial  spherical  shape  transforms  into  a fiat  superquadric.  Finally,  in 
figure  6.12(f),  all  the  points  on  our  shape  model  are  slopped,  thereby  recovering  the 
entire  shape  of  the  flat  superquadric.  Calculations  were  done  on  a 32  X 32  X 32  with 


Figure  6.4.  An  unsuccessful  attempt  to  reconstruct  a complex  shape  with  "signifi- 
cant" protrusions  usiug  an  active  contour  model.  Three  different  (poor)  results  are 
shown  in  parts  (c),  (f),  k (i)  corresponding  to  three  distinct  initialisations  in  parts 
(a),  (d),  k (g)  respectively. 


(c)  After  275  iterations  (f)  Alter  391  iterations 


Figure  B.5.  Reconstruction  of  a shape  with  ‘significant"  protrusions:  an  arterial  tree 
structure.  Computation  was  done  on  a 64  x 64  grid  with  a lime  step  At  = 0.00] . A 
narrow-band  updation  scheme  was  used  with  a band  width  of  6 = 0.075, 


(e)  Alter  125  iterations 


(1)  After  MO  iterations 


Figure  6.7.  Topological  split:  a single  instance  of  the  shape  model  splits  into  three 
instances  to  reconstruct  the  individual  shapes.  Compulation  was  done  on  a 64  x 64 
mesh  with  a time  step  At  = 0.001. 


Figure  6.8.  Topological  split  example  level  sets  shown  are  {0  = ±iC),  i € 1-5  +5| 
with  C = 0.(15. 
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(f)  Inner  boundaries 

Figure  6.9.  Shapes  with  holes:  a two-stage  scheme  is  used  to  arrive  at  a complete 
shape  description  of  both  simple  shapes  and  shapes  with  holes.  Computation  was 
performed  on  I2S  x I2S  grid  and  the  time  step  A I was  set  to  0.00025. 


(h)  Alter  10  iterations 


(c)  Alter  30  iterations 


Figure  6.10.  Expanding  torus  (F  = 1):  change  of  topology  is  depicted  in  this  example. 
Calculations  were  done  on  a 32  x 32  x 32  grid  with  a time  step  A I = 0.0005. 
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CONCLUDING  REMARKS 


while  retaining  the  desirable  features  of  existing  methods  for  shape  modeling,  over- 
comes most  of  their  deficiencies-  We  adopt  the  level  set  techniques  first  introduced 

recovery  which  rely  on  energy  minimisation,  the  final  result  in  our  method  is  com- 

featurcs  of  our  shape  modeling  scheme  include  its  ability  to  split  and  merge  freely 

tendibilily  to  higher  dimensions.  The  equations  of  motion  governing  our  evolutionary 
system  resemble  an  initial  valued  Hamillon-Jacobi  equation  with  a parabolic  right- 

Thus,  the  result  is  a very  genera)  shape  modeling  algorithm  which  we  believe  will 


In  the  context  of  2D  shape  recovery,  we  force  our  front  to  come  to  a stop  in 
images.  An  extension  is  defined  to  attribute  a consistent  meaning  to  this  speed 


complexity  is  realised  by  adopting  the  narrow-band  update  strategy.  In  the  narrow- 
band  scheme,  the  front  is  advanced  by  updating  the  0 function  in  a narrow  band 
around  the  zero  set.  The  efficacy  and  versatility  of  our  method  is  demonstrated  by  a 
series  of  shape  recovery  experiments  on  a set  of  synthetic  and  noisy  medical  images. 

Our  future  research  efforts  will  be  focussed  on  extending  this  method  to  other 
application  domains.  We  have  implemented  a similar  scheme  to  reconstruct  the 
surface  structures  of  3D  objects  from  3D  medical  image  data.  However,  it  is  not 
clear  how  to  adapt  the  level  set  formulation  iu  its  present  form  to  reconstruct  surface 
shapes  from  sparse  range  data.  The  issue  of  time  complexity  is  an  important  one, 
more  so  in  3D.  We  have  seen  in  chapter  4 that  the  stability  requirement  for  solving  the 
level  set  equations  forces  a very  small  time  step  for  line  grids  and  that  the  evaluation 


of  the  extension  is  computationally  very  expensive.  The  narrow-band  strategy  has 
drastically  reduced  the  computational  cost  in  2D,  but  implementing  a similar  scheme 
in  3D  could  prove  cumbersome.  One  approach  is  to  embed  our  level  set  algorithm 
in  a multiresolution  framework,  where  the  loss  of  accuracy  at  rapidly-converging 
coarse  grid  computation  is  compensated  by  highly  accurate  slowly-converging  fine 
grid  calculation.  Alternately,  sophisticated  parallel  implementations  can  be  devised. 
One  such  scheme  is  briefly  discussed  in  chapter  4. 
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